Isometric embeddings of black hole horizons in three-dimensional flat space 
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The geometry of a two-dimensional surface in a curved space can be most easily visualized by 
using an isometric embedding in flat three-dimensional space. Here we present a new method for 
embedding surfaces with spherical topology in flat space when such a embedding exists. Our method 
is based on expanding the surface in spherical harmonics and minimizing for the differences between 
the metric on the original surface and the metric on the trial surface in the space of the expansion 
coefficients. We have applied this method to study the geometry of back hole horizons in the presence 
of strong, non-axisymmetric, gravitational waves (Brill waves). We have noticed that, in many cases, 
although the metric of the horizon seems to have large deviations from axisymmetry, the intrinsic 
geometry of the horizon is almost axisymmetric. The origin of the large apparent non-axisymmetry 
of the metric is the deformation of the coordinate system in which the metric was computed. 
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I. INTRODUCTION 

In the past few years, important progress has been 
made in the development of fully three-dimensional (3D) 
numerical relativity codes. These codes are capable of 
simulating the evolution of strongly gravitating systems, 
such as colliding black holes and neutron stars, and can 
provide important physical information about those sys- 
tems such as the gravitational waves they produce. The 
codes also allow one to locate and track the evolution of 
apparent and event horizons of black holes that might 
exist already in the initial data or might form during 
the evolution of the spacetime. However, since the lo- 
cation of such horizons is obtained only in coordinate 
space, one typically has little information about the real 
geometry of those surfaces. One can, for example, obtain 
very similar shapes in coordinate space for horizons that 
are in fact very different (see for example the family of 
distorted black holes studied in [1] ; their coordinate loca- 
tions are very similar, but their geometries are quite dif- 
ferent). The most natural way to visualize the geometry 
of a black hole horizon, or of any other surface computed 
in some abstract curved space is to find a surface in or- 
dinary flat space that has the same intrinsic geometry 
as the original surface. The procedure of finding such a 
surface is called embedding the surface in flat space. 

It is a well known fact [2] that any two dimensional 
surface is locally embeddable in flat 3D space. Global 
embeddings of a surface, on the other hand, might easily 
not exist, and even when they do they are not easy to 
find. However, several methods have been proposed in 
the past for computing partial (or global) embeddings of 
surfaces when such embeddings exist. 

Partial embeddings of a slice through the Misner initial 



data for colliding black holes [3] have been computed for 
example in [4] . The method used to find such embeddings 
starts from the metric of the original surface written in 
terms of some local coordinates (u,v): 



ds 2 = Edu 2 + 2Fdudv + Gdv 2 



(1) 



One then introduces embedding functions X(u,v), 
Y(u,v) and Z(u,v) such that 

dX 2 + dY 2 + dZ 2 = ds 2 = Edu 2 + 2Fdudv + Gdv 2 , 



which implies 

E = 
F = 
G = 
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The above system of non-linear first order partial dif- 
ferential equations is not of any standard type. In order 
to solve it one can use a method originally proposed by 
Darboux. This leads to a single non-linear second or- 
der partial differential equation for Z[u, v) known as the 
Darboux equation. The character of the Darboux equa- 
tion depends on the sign of the Gaussian curvature K 
of the surface and on the orientation of the embedding 
at the point of integration as follows: for K > the 
equation is elliptic, for K < it is hyperbolic and it has 
parabolic character both if K = or if the surface is ver- 
tical. For the Misner geometry, the curvature is always 
negative and the Darboux equation is hyperbolic. It can 
then be rewritten by using the characteristics as coordi- 
nates, and solved as a Cauchy problem given appropriate 
initial data. Once the Darboux equation has been solved 
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for Z(u, v), the remaining equations can be integrated to 
give X(u,v) and Y(u,v). 

This method is not appropriate for embedding surfaces 
that have both regions with K > and regions with 
K < 0. The reason for this is that the Darboux equation 
one has to integrate is elliptic for K > and hyperbolic 
for K < 0. This means that for such surfaces the Dar- 
boux equation will change type and its integration will 
become very difficult. This will typically be the case for 
surfaces with spherical topology, such as black hole hori- 
zons, which will always have some regions of positive cur- 
vature, and may well have regions of negative curvature 
too. 

One of the first studies of the intrinsic geometry of 
rotating black hole horizon surfaces was carried out by 
Smarr [5]. There it was show, by direct construction of 
the embedding from the analytic Kerr metric, that while 
the horizon of a Schwarzschild black hole is spherical, for 
rotating black holes the horizon has an equatorial bulge, 
a satisfying and intuitive result that reinforces the no- 
tion that geometric studies of black hole horizons can add 
physical insight. The equatorial bulge can be character- 
ized by an oblateness factor that is uniquely determined 
by the ratio a/m, where m is the mass of the black hole 
and a its rotation parameter. It was also shown that 
for rapidly rotating black holes, with a/m > V3/2, the 
Gaussian curvature becomes negative near the poles, and 
the surface is not embeddablc in Euclidean space, as it is 
"too flat" ! 

Extending on this work, using a direct constructive 
embedding method described below, a number of studies 
were made of distorted, rotating, and colliding black hole 
horizons in axisymmetry [6,1,7-10] where it was shown 
that embeddings are very useful tools to aid in the un- 
derstanding of the dynamics of black holes. For exam- 
ple, distorted rotating black hole horizons were found to 
oscillate, about their oblate equilibrium shape, at their 
quasi-normal frequency. The recent work on isolated 
horizons [11-14] shows how geometric measurements of 
the horizon can be used to determine, for example, the 
spin of a black hole formed in some process, and other 
physical features. 

However, in the absence of axisymmetry, the problem 
of constructing an embedding for a black hole horizon 
becomes much more difficult. One approach to compute 
such embeddings of horizons in 3D spacetimes has been 
suggested by H.-P. Nollert and H. Herold [15]. They con- 
sider a triangular wire frame on the original surface and 
compute the distances between each point and its neigh- 
bors using the intrinsic metric of the surface. They then 
consider a network with the same topology in flat space 
and try to solve the system of equations 

\r i -r j \=d ij , (6) 

where represents the position vector of the k th point 
in flat space and dij represents the distance between the 



points Pi and Pj computed on the original surface. If 
necessary, they refine the grid until they reach a desired 
accuracy. 

The approach of Nollert and Herold seems very nat- 
ural, but is has the serious drawback that it does not 
always converge to the correct solution. The reason for 
this is that the method imposes constraints only on the 
distances between points, but it does not guarantee that 
the final surface will be smooth. There are in fact mul- 
tiple solutions to the system of equations, and for most 
such solutions the resulting embedding is not smooth. 
For example, if one tries to embed a simple sphere, this 
method might indeed converge to the sphere, but it might 
also converge to the surface one obtains if one cuts the 
top of the sphere, turns its up-side down and glues it 
back. The distances between point are the same in both 
surfaces, but only one of them is smooth. 

The method for computing embeddings that we 
present in this paper is based on a spectral decomposition 
of the surface in spherical harmonics written in a non- 
trivial mapping of the coordinate system. We search for 
the embedding by minimizing for the difference between 
the metric of the original surface and that of our trial em- 
bedding in the space of the coefficients of the spherical 
harmonics and of the coordinate mappings. Since we use 
a decomposition of the surface in spherical harmonics, 
the surface is guaranteed to be smooth. By increasing 
the number of spherical harmonics used in the decompo- 
sition of the physical surface, one can get as close to the 
correct embedding as desired. 

II. OUR METHOD 

The intrinsic geometry of any surface is completely de- 
termined by its metric. To construct the embedding of 
a given surface S, one needs to find a surface S' in flat 
space that has the same metric as S in an appropriate 
coordinate system. It is important to stress here the fact 
that finding the embedding surface S' also requires that 
one finds an appropriate mapping of the original coordi- 
nate system in the surface S to a new coordinate system 
in the surface S' in which the two metrics are supposed 
to agree. 

A. A Direct Method for Horizon Embeddings in 
Axisymmetry 

Before describing our new method, it is instructive to 
describe a simple and direct method for axisymmetric 
embeddings of horizons, used by a number of authors to 
study the physics of dynamic black holes [16,1,7-10,17]. 

For the case of non-rotating, axisymmetric spacetimes 
(easily generalized to rotation, but restricted here merely 



2 



for ease of illustration), the 3D metric on a given time 
slice can be written as 



gij{9, 4>) = ei{6, 4>) ■ e,-(0, <j>) = h ab e? e) , 



(10) 



A(r], 6)dr] 2 + B(r], 9)d9 2 + D(rj, 9) sin 2 9d<j) 2 . (7) 



The location of the axisymmetric horizon surface is given 
by the function n = rj s (9) . The 2D metric induced on the 
horizon surface is then given by 



ds 2 2 ) 



B + 



drjs 
d9 



A 



d9 2 + D sin 2 9d(j) 2 . (8) 



Now, the flat metric in cylindrical coordinates (z, p, ip) 
can be written as 



ds 2 = dz 2 + dp 2 + p 2 d^\f 



(9) 



To create an embedding in a 3D Euclidean space, we want 
to construct functions z(9, </>), p(9, (j)), and ip(9,(j)) such 
that we can identify the line elements given by Eqs. (8) 
and (9), that is, that all lengths be preserved. 

Here, we are faced with our first choice about the co- 
ordinates used in the embedding, a problem which will 
be more complex in the general case as we show below. 
Since the spacetime itself is axisymmetric, it is a nat- 
ural choice to make the embedding axisymmetric. We 
choose, then, to construct a surface for a constant value 
of <fi = 0, and we then have z — z{9), p = p(9), and the 
resulting embedding will be a surface of revolution about 
the z— axis. Using the obvious mapping between tp and 
(j), tp = <p, it becomes straightforward to derive ordinary 
differential equations to integrate for p{9) and z{9) along 
the horizon surface. It is important to emphasize that we 
have to make a choice about the embedding coordinates, 
even in this simpler case, as we must in the general case 
discussed below. 

Using this method, embeddings were carried out dur- 
ing the numerical evolution for a variety of dynamic, 
axisymmetric black hole spacetimes [16,1,7-10,17]. The 
evolution of these embeddings were found to be extremely 
useful in understanding the physics of these systems. We 
will use some of these results as test cases for the more 
general method for 3D spacetimes, as detailed in the next 
section. 



B. Our General Method for Embeddings in full 3D 

Given a coordinate system (i,j — 1,2) on our sur- 
face, the first step in looking for an embedding is to 
find the two dimensional metric gij of the surface in- 
duced by the metric of the three dimensional space h a b 
(a, b = 1,2, 3) in which it is defined. The general proce- 
dure to find such induced metric is to construct a coor- 
dinate basis of tangent vectors ei := di on the surface. 
The induced metric will then be given by 



where ef is the component of the vector with respect 
to the three dimensional coordinate x a . 

Since the surfaces we are concerned with in this paper 
(black hole horizons) have spherical topology, we will as- 
sume that the three-dimensional metric h a b is given in 
terms of spherical coordinates (r, 9, (f>) defined with re- 
spect to some origin enclosed by the surface. Further- 
more, we will also assume that the surface is a "ray- 
body" (Minkowski's strahlkorper [18], also known as a 
"star-shaped" region) that is, a surface such that any 
ray coming from the origin intersects the surface at only 
one point. Such a property implies that we can choose 
as a natural coordinate system on the surface simply the 
angular coordinates (9,(j)). 

If we take our surface to be defined by the function 
r = f(9,(f>), then it is not difficult to show that the in- 
duced metric on the surface will be given in terms of 
the three-dimensional metric h a t as: 

gee = hgg + h rr (dgff + 2h r8 dgf , (11) 
500 = ^00 + h rr (dtpf) 2 + 2h r< j,d^f , (12) 
ge<t> = hg^ + h rr dgfd^f + h r( pdgf + Kgd^f . (13) 

Let us now for moment a assume that an embedding of 
our surface in flat space exists, and let us also introduce a 
spherical coordinate system (r e , 9 e ,(j) e ) in flat space. No- 
tice that there is no reason to assume that a point with 
coordinates (9, (j)) in the original surface will be mapped 
to a point with the same angular coordinates in the em- 
bedding. In general, the embedded surface in flat space 
will be described by the relations 

r e =r e (0,0), 9 e = 9 e {9,4>), <j> e = <j> e {9, <f>) . (14) 

A crucial observation at this point is that the angu- 
lar coordinates {9, <f)} in the original surface still provide 
us with a well behaved coordinate system in the embed- 
ded surface, only one that does not correspond directly 
to the flat space angular coordinates {9 e ,<f) e }, but is in- 
stead related to them by the coordinate transformations 
8 e = 9 e (9,(f>) and (f) e — (f> e {9,4>). This means that un- 
der the embedding, points with coordinates (9, </>) in the 
original surface will be mapped to points with the same 
coordinates (9, </>) in the embedded surface, but different 
coordinates (9 e ,(j) e ). We then have two natural sets of 
coordinates in the embedded surface: the ones inherited 
directly from the original surface through the embedding 
mapping, and the standard angular coordinates in flat 
space. 

By definition, an embedding preserves distances, so the 
proper distance between two points in the original sur- 
face must be equal to the distance between the two cor- 
responding points in the embedded surface. Since those 
corresponding points have precisely the same coordinates 
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{9, </>}, we must conclude that for the embedding to be 
correct, the components of the metric tensor in both sur- 
faces when expressed in terms of the coordinates {9, </>} 
must be identical. That is, if we call gfj the metric of the 
embedded surface, we must have 



gee = . 



9e<t> — ge<p > 9<t>4> ~ 9<t><t> ■ 



(15) 



It is important to stress that the components of the met- 
ric in the embedded surface will only be equal to the 
components of the metric in the original surface if we use 
the inherited coordinate system, but not if we use the 
standard angular coordinates in flat space. 

Computing the components of the metric in the em- 
bedded surface in terms of the inherited coordinate sys- 
tem {9,(j)}, given the embedding relations (14), is not 
difficult. All one needs to do in practice is consider four 
points in the original surface with coordinates P\{9,<j>), 
P 2 = (9 + 59, <f>), P 3 = (6,<f> + 50), Pi = (9 + 59, 4> + 8<j>), 
find their corresponding coordinates {r e ,9 e ,(f> e } in flat 
space using (14), compute their squared distances using 
the flat space metric, and then solve for the metric com- 
ponents from 



(PiP 3 f = g^dtf , 

(Pj\) 2 = gle d9 2 + gfo d<f? + 2g* H d9d0 . 



(16) 
(17) 
(18) 



Finding the embedding now means finding a map- 
ping (14) such that equations (15) are satisfied every- 
where. 

Let us consider first the relation between the (9, <f>) 
coordinates on the original surface and the angular coor- 
dinates (9 e , <j> e ) in flat space. Even if these two sets of co- 
ordinates are not equal, we can safely assume that there 
is a one-to-one correspondence between them. Moreover, 
both are sets of angular coordinates, so they have the 
same behavior: 9 and 9 e go from to w, and <p and <j) e 
go from to 2it and are periodic. From these properties, 
it is not difficult to see that the most general functional 
relation between both sets has the form 

oo 

9 e {0,4) = 9 + J2bnosm(n9) 

oo oo 

+ ^2 ^2 b nm sm(nff) sin(m0) , (19) 

n—l m—1 

oo oo 

<p e (9,(f>) = <f> + cqo + ^ c„ sin(n#) + ^ c 0m sin m<f) 

n—l m—1 

oo oo 

+ c nm sm(n9) sin(m0) . (20) 



n—l m—1 



The second term in the expansion for 9 e represents a 
general axisymmetric re-mapping of 9, while the third 
term is required if axisymmetry is not assumed. In the 



expression for <f> e , the second and third terms represent 
a possible rigid twist of the coordinate system, and the 
last two terms stand for a general dependence of <p e on 
both 9 and <ft. 

For the radial coordinate r e , it is also not difficult to 
see that one can use a simple expansion in spherical har- 
monics of the form 



^ ^ V4^ai m Yi m (9,(f)) 



(21) 



Z=0 m=-l 



where the overall normalization factor of v47r has been 
inserted so that aoo is the average radius of the surface, 
aio is its average displacement in the z-direction, and so 
on. We will also use a real basis of spherical harmonics, 
for which m and — m stand for an angular dependence 
cos(to</>) and sin(m0), instead of the complex exp(im^) 
and exp(— im<p). 

The metric of the embedded surface will now be com- 
pletely determined by the set of coefficients ai m , b nm and 
c nm . The space of these coefficients can be regarded as 
a vector space V, with any given point in V representing 
a surface in flat space together with a certain coordinate 
mapping. 

Consider now an embeddable surface S in some arbi- 
trary curved space. It is not difficult to find a real valued 
function F defined on V that has a global minimum at 
the point P € V for which the metric of the embedded 
surface S e is the same as the metric of the original sur- 
face. One such function is 



/>7T />2"7T 

J 8=0 J(b=0 



{gee{9, 



>(M)r 



+ {g^{9,4>)~g e ^{9,^))' 



+ 



4>)-9'h(O,0)) 2 d9dcf>. (22) 



We call the function F above the "embedding" function. 
It is easy to see that F > on any point in V, and that 
F = if and only if g — g e for all (9, 4>). The embedding 
then corresponds to the absolute minimum of F in V. 
There are many different numerical algorithms for finding 
minima of general functions in multidimensional spaces. 
In our code we have used Powell's minimization algo- 
rithm [19], but we are aware that other methods might 
perform better. It is important to mention that the defi- 
nition of the embedding function F above is by no means 
unique. Many different forms for F can be constructed, 
in particular, one could take into account the fact that 
not all metric functions have similar magnitudes and con- 
struct an embedding function that normalizes each term 
in the above expression. 

One of the problems with minimization algorithms in 
general is that they cannot distinguish between a global 
minimum (what we really want) and a local minimum. 
In our case, the value of F at the absolute minimum 
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is zero, so we can easily distinguish between a real em- 
bedding and a wrong solution that might appear if the 
minimization algorithm gets stuck in a local minimum. 
However, steering the algorithm toward the global min- 
imum is non-trivial. From experience we have seen that 
local minima for which F ^ do exist for our problem. 
In order to avoid them, we have found it necessary to run 
first the minimization algorithm with a small number of 
coefficients a, b, and c, and then increase the number of 
coefficients one by one until we find a good solution. This 
method is certainly time consuming, but it seems to work 
well in the examples we have considered so far. 

Of course, in order to find a "perfect" embedding, one 
would have to push the number of coefficients all the way 
to infinity. This is numerically impossible, so in practice 
we just set up a given tolerance in the value of the func- 
tion F and increase the number of coefficients until we 
achieve that tolerance. We also check that the value of F 
goes to zero exponentially as we increase the total num- 
ber of coefficients. We have seen that n = I <~ 14 is 
enough for relatively simple surfaces like most black hole 
horizons. If one wants to embed something more com- 
plicated (like a human face, for example) starting from 
its metric in some coordinate system, this value would 
clearly be too small. Another way to see whether an em- 
bedding is good or not is to compare directly the metric 
of the original surface with the metric of the resulting 
embedding. If the fit is good enough, the code has con- 
verged to the correct embedding. 

One important test we have used for our algorithm is 
a direct comparison of the results obtained with our code 
with embeddings computed with a different code in the 
special case when the surface is axisymmetric. 
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2 


2 + 1.9 x 1(T 5 



TABLE I. Comparison of the recovered expansion coeffi- 
cients for the embedding of the test surface described in the 
text. 

the lower panel the resulting embedding. The differences 
in the shape of the two surfaces are very difficult to see. 

For this test we have used 100 x 100 grid points to de- 
scribe the surface. Since the surface is symmetric with re- 
spect to reflections on all three (x, y, z) coordinate planes, 
we have considered only one octant, so the angular reso- 
lution was A9 — Ac/) = 7r/200. The recovered expansion 
coefficients are shown in Table I. The coefficients cor- 
responding to the mapping of the angular coordinates, 
as well as the rest of the Yi m coefficients were either ex- 
actly zero because of the octant symmetry or had values 
smaller than 10~ 3 . 

Figure 2 shows a direct comparison of the metric com- 
ponents ggg , and along the lines 9 = 7r/4 and 
4> = 7r/4, i.c in the middle of the computational domain. 

The final value for the embedding function in this test 
was F = 1.6 x 10~ 5 , but we have found that we can easily 
decrease this value by refining the numerical grid on the 
surface. 



III. TESTS 
A. Recovering a known surface 

A very simple test for our algorithm is to look for the 
embedding of a surface that is known to be embeddable 
and has a known embedding. In order to do this we first 
construct a surface in flat space by choosing some arbi- 
trary values of the spherical harmonic coefficients. We 
then compute the metric of this surface in the standard 
(9, (f>) angular coordinates, and give this metric as input 
to our code. The code must then recover the correct val- 
ues of the spherical harmonic coefficients plus a trivial 
mapping of the angular coordinates. 

We show an example of this in Fig. f , where we have 
chosen a surface defined by the spherical harmonic coef- 
ficients: 

«oo = 9 , a 22 = 1 , a 44 = 4 , (23) 

with all other coefficients equal to zero. In the upper 
panel of the figure we show the original surface, and in 



B. An axisymmetric example: Rotating black holes 

As already mentioned in section I, a well-known set 
of axisymmetric surfaces whose embeddings have been 
studied, first by Smarr [5], and also as a test case in [15], 
is that of the horizons of rotating black holes. In the 
static Kerr case, the metric of the horizon is given by 

da 2 = p 2 d9 2 + ^ (r 2 + a 2 ) d<t> 2 , (24) 

with 

p 2 = r 2 + a 2 cos 2 9 , r = m + \/ m 2 — a 2 , (25) 

and where m and a are two parameters representing the 
mass of the black hole and its angular momentum respec- 
tively. 

It is well known that the Kerr horizon is globally em- 
beddable in flat space only for a/m < [5]. Using 
our embedding code, we have been able to successfully 
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FIG. 1. Embedding of a test surface defined by the spher- 
ical harmonic coefficients (aoo = 9, 022 = 1,144 = 4). The 
upper panel shows the original surface and the lower panel 
the resulting embedding. 



recover these embcddings when they exist. As an exam- 
ple we show in Fig. 3 the embedding obtained for the last 
embeddable case a/m = y/3/2. The solid line shows the 
embedding obtained with the axisymmetric algorithm de- 
scribed above, and the dotted line the one obtained with 
our minimization algorithm. 

When a/m > Vo/2, a global embedding does not ex- 
ist, and our method fails as expected. Figure 4 shows 
the result of an attempt to embed a Kerr black hole hori- 
zon in the case when a/m — 0.99. The dotted line is 
the output of our minimization algorithm and the solid 
line is an embedding of the same surface made with an 
axisymmetric algorithm in the embeddable region, plus a 
flat top in the region where the embedding does not exist. 
The axisymmetric method is a local constructive method, 
and hence it is able to build the embedding surface from 
one point to the next where it exists (in this case start- 
ing from the equator). Since our method is global we get 
the embedding wrong everywhere. As currently imple- 
mented, our method will insist on trying to find a global 
embedding, and will settle on a shape that minimizes the 
function F. If the embedding does not exist, the mini- 
mum value of F found will be clearly different from zero. 
This is easily seen by examining the residual function F 
for the embeddings of the Kerr horizons. In Fig. 5 we 
show the value of the minimum of F found with our al- 
gorithm, plotted against the total number of expansion 
coefficients, for the cases a/m = \/3/2 (solid line) and 
a/m — 0.99 (dotted line). One can see that for the non- 
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Minimization algorithm 

Axisymmetric algorithm — embeddable region 



FIG. 3. Embedding of a Kerr black hole with a/m — v3/2. 
The dotted line is the embedding we obtained using our min- 
imization algorithm. The solid line is an embedding of the 
same surface computed with an axisymmetric algorithm. 



embeddable case, F stops decreasing at a value that is 
more than 2 orders of magnitude larger than the one we 
obtain when the embedding exists. 

Our code might be adjusted in the future for finding 
partial embeddings by reducing the integration domain 
in the definition of F, Eq. (22), to something smaller 
than the full sphere. This approach might lead to correct 
partial embeddings of surfaces that cannot be embedded 
globally. The disadvantage would be that one would have 
to guess the domain where the embedding exits before 
starting the computation. 



C. Black hole plus Brill wave 

We now move to the case of numerically generated, dis- 
torted black holes. Schwarzschild black holes distorted 
by Brill waves [20] have been extensively studied in nu- 
merical relativity [21-24]. The axisymmetric data sets 
used for numerical evolutions consist of a Schwarzschild 
black hole distorted by a toroidal, time-symmetric grav- 
itational wave. 

It is convenient to describe the metric of the black 
hole plus Brill wave spacetime in a spherical-polar like 
coordinate system (77, 9, <f>) were 77 is a logarithmic radial 
coordinate defined by rj — ln(2r/M) and (0,(f>) are the 
standard angular coordinates. In these coordinates, the 
spatial metric has the form [6,23] 



FIG. 4. An attempt to embed a Kerr black hole horizon 
with a/m — 0.99. The dotted line is the output of the min- 
imization algorithm. The solid line is an embedding of the 
same surface made with an axisymmetric algorithm. The flat 
line on the top represents the region where an embedding in 
flat space does not exist. 




v-v 



last embeddable Kerr horizon 



non— embeddable Kerr with a/m=0.99 



10 100 
total number of coefficients 



FIG. 5. The value of the embedding function F versus the 
total number of coefficients for the two Kerr black holes dis- 
cussed above. The triangles correspond to a/m = 0.99 and 
the stars to a/m — V3/2. 
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dl 2 = * 4 [e 2 " (d V 2 + d9 2 ) + r 2 sin 2 Odcf 



(26) 



where both q and are functions of rj and only. In or- 
der to satisfy the appropriate regularity and fall-off con- 
ditions, the function q has been chosen in the following 

way 



q{rj,6) = a sin™ (9 



e { — > + e 



(27) 



with n is an arbitrary even number larger than 0. The pa- 
rameter a characterizes the amplitude of the Brill wave, 
while the parameters b and to characterize its radial loca- 
tion and width respectively. Having chosen the form of 
the function q, the hamiltonian constraint is solved nu- 
merically for the conformal factor ty. An isometry condi- 
tion is imposed at a coordinate sphere to guarantee that 
the final spacetime will contain a black hole. 

Notice that the metric (26) is the three-dimensional 
metric of space, and not the two-dimensional metric of 
the apparent horizons. The apparent horizons for these 
data sets have to be located numerically. Once these 
horizons are found, their two-dimensional metric can be 
computed from the three-dimensional metric given above, 
using the expressions given in (13). 

The horizons from these axisymmctric black hole plus 
Brill wave data sets and their embeddings have been 
studied previously in [23] and we have been able to re- 
produce their results using our algorithm. An example of 
this can be seen in Fig. 6, where we show the embedding 
of the horizon of a black hole plus Brill wave data set 
corresponding to the parameters 



a =1.0, 6 = 0.0, 



= 1.0. 



(28) 



In the figure, the dotted line shows the embedding ob- 
tained with our minimization algorithm, and the solid 
line shows the embedding of the same surface obtained 
in [23] . Notice how the intrinsic geometry of the horizon 
is far from spherical. 



D. Application to Full 3D Spacetimes 

Having tested our algorithm on both analytic and nu- 
merically generated axisymmetric black hole spacetimes, 
we now turn to the case of full 3D black hole spacetimes 
for which our method was developed. 

The axisymmetric black hole plus Brill wave data sets 
from [23] have been generalized in [25] to full 3D by mul- 
tiplying the Brill wave function q by a factor that has 
azimuthal dependence to obtain 



q(r], 8,4>) = a sin™ 9 (l + c cos 2 (j>) 
-(^) 2 ) +e -(^) 21 



(29) 



were c is an arbitrary parameter characterizing the non- 
axisymmetry of the Brill wave. 




FIG. 6. Embedding of the apparent horizon of a black 
hole plus Brill wave data set corresponding to the param- 
eters (a = 1.0,6 = 0.0, lu = 1.0, n = 2). The dotted line is 
the embedding obtained with our minimization algorithm and 
the solid line the embedding of the same surface obtained by 
Anninos et al. 



We have computed embeddings of the non axisymmct- 
ric apparent horizons obtained in this case. Here we will 
show examples of two such horizons. First we consider 
the embedding of the apparent horizon for the data set 
with parameters 

a = 0.3, 6 = 0.0, u = 1.0, n = 4, c = 0.4, (30) 

which corresponds to a relatively small non-axisymmctric 
distortion of the black hole. Figure 7 shows the em- 
bedding of the corresponding apparent horizon. Notice 
how the surface looks quite axisymmetric even though 
we have added a non-trivial non-axisymmctric contribu- 
tion to the metric. It is clear that the non-axisymmctry 
of the metric components is to a large degree a coordi- 
nate effect. Numerical evolutions of such black holes do 
show radiation in non-axisymmetric modes of gravita- 
tional radiation. However the mass energy carried away 
by the non-axisymmetric modes is much smaller than the 
energy of the axisymmctric modes [26,21]. This is con- 
sistent with our result showing that the horizon is almost 
axi-symmetric. 

Figure 8 shows a direct comparison of the different 
angular metric components on the apparent horizon and 
the resulting embedding, along the </> = 7r/4 and 6 = 7r/4 
lines. We can see how the fit is very good in both cases. 
Notice also how there is indeed some dependence of the 
metric components on (p. 

To check if our algorithm is converging to the correct 
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FIG. 7. Embedding of the apparent horizon for the non ax- 
isymmetric black hole plus Brill wave data set corresponding 
to the parameters (a = 1.0,6 = 0.0, u) = 1.0, n = 4, c = 0.4). 
Although the metric has a non-trivial non-axisymmetric con- 
tribution, the surface looks quite axisymmetric. 



40 




_5 i ^ , ^ Lio i ^ ■ ^ 1 

it/4 jt/2 it/4 nil 

FIG. 8. We show the angular metric components on the 
apparent horizon and on the resulting embedding for the black 
hole plus Brill wave data set corresponding to the parameters 
(a = 1.0, b = 0.0, w = 1.0, n = 4,c = 0.4). On the left panel 
we show the line 9 = 7r/4 and on the right panel the line 

(j> = 7T/4. 



1e-02 



1e-03 r 



1e-04 r 



1e-05 r 




1e-06 1 ■ — ■ — ■ — 

1 10 100 1000 

total number of coefficients 

FIG. 9. We show the value of the embedding function 
F at the minimum in terms of the total number of expan- 
sion coefficients in a logarithmic scale for the black hole 
plus Brill wave data set corresponding to the parameters 
(a = 1.0, 6 = 0.0, w = 1.0, n = 4, c = 0.4). 

embedding, we show in Fig. 9 the value of the embedding 
function F at the minimum, in terms of the total number 
of expansion coefficients. We clearly see that the value 
of F is converging exponentially to zero. 

As a second example, we now consider the embedding 
of the apparent horizon corresponding to the black hole 
plus Brill wave data set with parameters 

a = 0.3, 6 = 0.0, u = 1.0, n = 4, c=1.9. (31) 

In this case, the non-axisymmetry is considerably larger, 
and one can see from Fig. 10 that the horizon is clearly 
not axisymmetric. 

In Fig. 11 we show again a direct comparison of the 
angular metric components on the apparent horizon and 
the resulting embedding along the <j> = 7t/4 and = 7r/4 
lines. As before, the fit is very good. 

Finally, in Fig. 12 we show again the value of the em- 
bedding function F at the minimum in terms of the num- 
ber of expansion coefficients. As before, the value of F 
converges exponentially to zero, but the convergence is 
slower than in the previous example due to the higher 
degree of complexity of the surface. 

IV. CONCLUSIONS 

We have implemented and tested a new algorithm for 
computing isometric embeddings of curved surfaces with 
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FIG. 11. The three independent components of the 
metric for the black hole plus Brill wave will a larger 
non-axisymmetric perturbation. On the left panel we show 
the line 6 = it/4 and on the right panel the line <f> — tv/4. 




u- 1 e-04 




FIG. 10. Two orientations of the embedding of the appar- 
ent horizon of a black hole perturbed by a Brill wave with a 
higher non-axisymmetry than in the previous example. 



10 100 
total number of coefficients 



FIG. 12. The value of the embedding function F at the 
minimum in terms of the total number of expansion coeffi- 
cients on a logarithmic scale for the black hole plus Brill wave 
with a larger non-axisymmetric contribution. 
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spherical topology in flat space. We define a function on 
the space of surfaces that has a global minimum for the 
right embedding and we find this minimum using a stan- 
dard minimization algorithm. In this paper we have dis- 
cussed the method and its applications to black hole vi- 
sualization in numerical relativity. The method has been 
tested both on simple test surfaces and on Kerr black 
hole horizons, and shown to correctly reproduce known 
results within specified tolerances. We have also used our 
method construct embeddings of non-axisymmctric, dis- 
torted black holes for the first time. We observed that 
the non-axisymmetry of the embedded surface is some- 
how smaller than one expects from just looking at the 
metric. This is consistent with the small amount of grav- 
itational radiation emitted in non-axisymmetric modes 
during the numerical evolution of such systems [21,26]. 

Our method is rather robust, and by construction pro- 
duces smooth surfaces, therefore avoiding some of the 
problems of previous methods. One disadvantage of our 
method is that the expansion in spherical harmonics im- 
plies that it can only be used to embed ray-body surfaces 
(i.e. surfaces such that any ray coming from the origin 
intersects the surface at only one point), still we expect 
most black hole horizons to have this property. The main 
disadvantage, however, is that the method is very time 
consuming due to the fact that minimization algorithms 
in general are slow. Another problem is the fact that 
minimization algorithms can easily get trapped in local 
minima. In order to avoid this we have found it necessary 
to increase the number of coefficients one by one and to 
use at each step the result of the previous step as initial 
guess, adding to the total amount of time the algorithm 
needs to find the embedding. As presently implemented, 
it also cannot find partial embeddings when no global 
embeddings exist, but straightforward modifications to 
the algorithm should permit this in certain cases. 

In the future, we will apply this method to study the 
dynamics of 3D black hole horizons as a tool to aid in 
understanding the physics of such systems. Although our 
method has been applied in this paper only to apparent 
horizons, it can clearly be applied to obtain embeddings 
of event horizons as well, once they have been located in 
numerical evolutions. 

Our embedding algorithm has been implemented as a 
thorn in the Cactus code [27], and is available for the 
community upon request from the authors. 
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